###CCA
setwd("D:/入侵/数据/分析2021/5 CCA")
sp=read.table("otu.csv",header=T,sep=",",row.names=1)
sp_dat=sp[-1,]
s.env=read.table("env.csv",header=T,sep=",",row.names=1)
group=read.table("names.csv",header=T,sep=",",colClasses=c("character","character")) 

#cca
library(vegan)
cca=cca(sp_dat~.,s.env,scale=T) 
sam=data.frame(cca$CCA$u[,c(1,2)],group$Community) 
colnames(sam)=c("CCA1","CCA2","group")
sam$group = factor(sam$group, levels=c('OP','IP','IPS','IS','RP'))

env=cca$CCA$biplot[,c(1,2)]  %>% as.data.frame()
cca1=round(cca$CCA$eig[1]/sum(cca$CCA$eig)*100,2) 
cca2=round(cca$CCA$eig[2]/sum(cca$CCA$eig)*100,2) 

p=ggplot(data=sam,aes(CCA1,CCA2))+
  geom_point(aes(colour=group,shape=group),size=5)+
  scale_color_manual(values=c("#44BBBB", "#8ea1ca", "#E6BD1A", "#8CBB44", "#f28a6a"))+
  geom_text(aes(label=rownames(sam), group = group, colour = group),size=4,hjust=0.5,vjust=-0.7,position="jitter")+ 
  scale_shape_manual(values=c(19,19,19,19,19))+
  labs(x=paste("CCA1 ",cca1," %"),y=paste("CCA2 ",cca2," %"))+
  theme(text=element_text(family="serif"))+
  theme_customize+
  stat_ellipse(aes(group = group, colour = group),level = 0.95, show.legend = F)+
  geom_segment(data = env,aes(x=0,y=0,xend = env[,1]*1.8, yend = env[,2]*1.8), colour="#C43C3C", size=1, arrow=arrow(angle=35, length=unit(0.3, "cm"))) +
  geom_text(data=env, aes(x=env[,1]*1.8, y=env[,2]*1.8, label=rownames(env)),
            size=5, colour="#C43C3C",hjust = (1 - 2 * sign(env[ ,1])) / 3,
            angle = (180/pi) * atan(env[ ,2]/env[ ,1]))+ # #8*7
  geom_hline(yintercept=0, linetype=5) +geom_vline(xintercept=0, linetype=5) 

sink("CCA_ANOVA_result.txt") 
print(anova.cca(cca,step=999))
print(anova.cca(cca,by="term",step=999))
sink()

#####CCA-TROPHIC
tro_dat=read.table("otu_tro.csv",header=T,sep=",",row.names=1)
s.env_t=read.table("env.csv",header=T,sep=",",row.names=1)
group_t=read.table("names.csv",header=T,sep=",",colClasses=c("character","character"))

#cca
cca_t=cca(tro_dat~.,s.env_t,scale=T) 
sam_t=data.frame(cca_t$CCA$u[,c(1,2)],group_t$Community) 
colnames(sam_t)=c("CCA1","CCA2","group")
sam_t$group = factor(sam_t$group, levels=c('OP','IP','IPS','IS','RP'))

env_t=cca_t$CCA$biplot[,c(1,2)]  %>% as.data.frame()
cca1_t=round(cca_t$CCA$eig[1]/sum(cca_t$CCA$eig)*100,2)
cca2_t=round(cca_t$CCA$eig[2]/sum(cca_t$CCA$eig)*100,2)

p_t=ggplot(data=sam_t,aes(CCA1,CCA2))+
  geom_point(aes(colour=group,shape=group),size=5)+
  scale_color_manual(values=c("#44BBBB", "#8ea1ca", "#E6BD1A", "#8CBB44", "#f28a6a"))+
  geom_text(aes(label=rownames(sam), group = group, colour = group),size=4,hjust=0.5,vjust=-0.7,position="jitter")+ 
  scale_shape_manual(values=c(19,19,19,19,19))+
  labs(x=paste("CCA1 ",cca1_t," %"),y=paste("CCA2 ",cca2_t," %"))+
  theme(text=element_text(family="serif"))+
  theme_customize+
  stat_ellipse(aes(group = group, colour = group),level = 0.95, show.legend = F)+
  geom_segment(data = env_t,aes(x=0,y=0,xend = env_t[,1]*1.8, yend = env_t[,2]*1.8), colour="#C43C3C", size=1, arrow=arrow(angle=35, length=unit(0.3, "cm"))) +
  geom_text(data=env_t, aes(x=env_t[,1]*1.8, y=env_t[,2]*1.8, label=rownames(env_t)),
            size=5, colour="#C43C3C",hjust = (1 - 2 * sign(env_t[ ,1])) / 3,
            angle = (180/pi) * atan(env_t[ ,2]/env_t[ ,1]))+ # #8*7
  geom_hline(yintercept=0, linetype=5) +geom_vline(xintercept=0, linetype=5) 

sink("CCA_ANOVA_result_tro.txt") 
print(anova.cca(cca,step=999))
print(anova.cca(cca,by="term",step=999))
sink()

multiplot(p, p_t, cols=2) 

##bar graph of environment variables
##variance analysis
library(vegan)
library("reshape2")
env=read.table("env0.csv",header=T,sep=",")
dat1<-melt(env,
           id.vars = c('Community','Site'),
           variable.name='ind',
           value.name='val')

library(dplyr)
list_name1 <-  unique(dat1$ind)
##normality    ########leaf c/n &soil ph √   ####>0.05 normality
sink("normal-test result.txt") 
for(i in list_name1){
  print(paste(i))
  tmp<-shapiro.test((dat1 %>% filter(ind == i))$val)
  print(tmp)
}
sink()

##homogeneity of variance    ########leaf c/n &soil ph √    ####>0.05 homogeneous
sink("leveneTest result.txt")
for(i in list_name1){
  print(paste(i))
  tmp<-dat1 %>% filter(ind == i) %>% leveneTest(val,Community,data=.)
  print(tmp)
}
sink()

dat_test=dat1 %>% filter(ind == i)
leveneTest(val,Community,data=dat_test)

library(agricolae)
sink("aov result.txt") ##leaf n
for(i in list_name1){
  print(paste(i))
  mod<-dat1 %>% filter(ind == i) %>% aov(val~Community,data=.)
  print(summary(mod))
  out<- LSD.test(mod,"Community",p.adj="none" )
  print(out$groups)
}
sink()

sink("Kruskal result.txt") 
for(i in list_name1){
  print(paste(i))
  mod<-dat1 %>% filter(ind == i) %>% with(.,kruskal(val,Community,group=TRUE, main=paste(i,sep = " ")))
  print(mod)
}
sink()

env_bar=read.csv("env_bar.csv",header=TRUE)
list_env <-  unique(env_bar$ind)
for(i in list_env){
  tmp<-env_bar %>% 
    filter(ind == i) %>% 
    dcast(., Community ~ type, value.var = "val") 
  tmp$Community = factor(tmp$Community, levels=c('OP','IP','IPS','IS','RP'))
  tmp$mean=as.numeric(tmp$mean)
  tmp$SE=as.numeric(tmp$SE)
  assign(paste("p_", i, sep = ""), ggplot(tmp,aes(x=Community,y=mean,fill=Community))+
           theme_customize+
           theme(legend.position="none")+
           geom_bar(position="dodge",stat="identity",width=0.5)+
           geom_text(aes(y=mean+3*SE,label=text),position=position_dodge(width = 0.9))+ #Add text to the bar graph
           xlab("Plant Community")+
           ylab(paste(i))+
           geom_errorbar(aes(ymax = mean+1.96*SE, ymin = mean-1.96*SE), 
                         position = position_dodge(0.9), width = 0.15))
}

multiplot(p_Aboveground_biomass, 
          p_Plant_density, 
          p_Leaf_N,
          p_Soil_salinity,
          cols=4)  ##12*3